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ABSTRACT 

Context. Cosmic rays passing through large astrospheres can be efficiently cooled inside these “cavities” in the interstellar medium. 
Moreover, the energy spectra of these energetic particles are already modulated in front of the astrospherical how shocks. 

Aims. We study the cosmic ray flux in and around A Cephei as an example for an astrosphere. The large-scale plasma flow is modeled 
hydrodynamically with radiative cooling. 

Methods. We studied the cosmic ray flux in a stellar wind cavity using a transport model based on stochastic differential equations. The 
required parameters, most importantly, the elements of the diffusion tensor, are based on the heliospheric parameters. The magnetic 
field required for the diffusion coefficients is calculated kinematically. We discuss the transport in an astrospheric scenario with 
varying parameters for the transport coefficients. 

Results. We show that large stellar wind cavities can act as sinks for the galactic cosmic ray flux and thus can give rise to small-scale 
anisotropies in the direction to the observer. 

Conclusions. Small-scale cosmic ray anisotropies can naturally be explained by the modulation of cosmic ray spectra in huge stellar 
wind cavities. 
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1. Introduction 


Recently, simulations of astrospheres around hot stars have 
gained new interest, see for example |Decin et al.| (|20T2|l, |Cox| 
|et al.| P012| l, I Arthur] ( |2012| l, [van M^le et al.| ( |2014| l. These 
authors modeled astrospheres using a (magneto-)hydrodynamic 
approach, either in ID or 2D. In this work, such astrospheric 
models are used for the first time to estimate the cosmic ray flux 
(CRF) through it. Because of the large spatial extent of O star as¬ 
trospheres (wind bubbles), these objects can efficiently cool the 
spectrum of galactic cosmic rays (GCR). Especially A Cephei is 
an interesting example, being the brightest runaway Of star in 
the sky (type 06If(n)p). We estimate the CRF at different ener¬ 
gies for A Cephei as an example of an O star astrosphere using 
stochastic differential equations (SDE) ( jStrauss et al. [2013| l to 
solve the GCR transport equation. Runaway O and B stars are 
common and part of a sizable population in the galaxy, a signifi¬ 


cant number of which show bow-shock nebulae (e.g., Huthoff & 
|Kaper| [2002| l. In Sect. we discuss the radiative cooling func¬ 
tions, while in Sect.j^we show the astrosphere model results. In 
Sect. |5]we estimate the CRF. 


2. Large-scale structure of astrospheres 

Winds around runaway stars, or in general, stars with a 
nonzero relative speed with respect to the surrounding interstel¬ 
lar medium (ISM), develop bullet-shaped astrospheres. 


The hydrodynamic large-scale structure is sketched in Fig.[T] 
the notation of which is described below. 

The hypersonic stellar wind (Mach numbers Ma » 1) un¬ 
dergoes a shock transition to subsonic velocities at the termina¬ 
tion shock (TS) in the inflow direction. Then a tangential dis¬ 
continuity, the astropause (AP), is formed between the ISM and 
the stellar wind, where the velocity normal to it vanishes: there 
is no mass transport through the AP. Other quantities such as 
the tangential velocity, temperature, and density are discontinu¬ 
ous, while the thermal pressure is the same on both sides. If the 
relative speed, or the interstellar wind speed as seen upwind in 
the rest frame of the star, is supersonic in the ISM, a bow shock 
(BS) exists. If the relative speed is subsonic, there will be no BS, 
see Table for the stellar parameters and for the stellar-centric 
model distances of the TS, AP, and BS. The region between the 
BS and AP is called outer astrosheath, the region between the 
AP and TS the inner astrosheath. The AP around the inflow di¬ 
rection at the stagnation line is sometimes called the nose, while 
the region beyond the downwind TS is called the astrotail. The 
latter can extend deep into the ISM. The region inside the TS is 
called the inner astrosphere. 

In the downwind direction, the termination shock forms a 
triple point, from which the Mach disk (MD) extends down to 
the stagnation line; this is the line through the stagnation point 
and the star. A tangential discontinuity (TD) emerges down into 
the tail. A reflected shock (not shown here) also extends from 
the triple point toward the TD. 
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text. 


This is the standard shape of an astrosphere using a single 
hydrodynamic fluid ( Baranov et al.[ 1971[ Pauls et aT] 1995] l. 
For the large dimensions of O-star astrospheres, cooling oper¬ 
ates inside the outer astrosheath, but usually not in the inner as- 
trosheath. Cooling is also present beyond the model boundary 
(see Table of the inner astrosphere, which leads to relative 
sizes of these regions different from that of pure hydrodynamic 
flow. 

The astrosphere is always bullet-shaped, which can be seen 
by the conservation of momentum: PswV^w + Psw - Pismv^sm 
P iSMi where p,v,P are the density, velocity, and the thermal 
pressure of the stellar wind (subscript sw) and the ISM (subscript 
ism)- The stellar wind pressure is usually neglible inside the TS. 
For hypersonic flows this holds true for the thermal ISM pres¬ 
sure in inflow direction, while for subsonic inflows the thermal 
pressure dominates. In the tail direction, only the thermal ISM 
presssure is present beyond the TS. This means that even if the 
inflow is subsonic, there is a total pressure asymmetry between 
upwind and downwind as long as the ISM velocity does not van¬ 
ish. Thus on long time and large spatial scales, a bullet-shaped 
astrosphere will develop. 


3. Cooling and heating 


As a result of the large dimensions of astrospheres and their outer 
astrosheath, the plasma can effectively be cooled by Coulomb 
collisions. In this process, an electron in an atom (molecule) will 
be excited, and after returning to a lower energy state, the re¬ 
emitted photon will carry away the energy ([Sutherland & Dopita 


|1993| and references therein). This can lead to an effective cool 
ing of the shocked ISM (ISM 2 in Fig.[T]l. 

Several cooling functions are discussed in the literature (e.g., 
Rosner et al.j |1978[ Mellema & Lundqvistj 2002 1 |Townsend| 


2009 Schure et al. 2009 1 Reitberger et al. 


2014 1 . The differ¬ 


ences in the cooling functions are caused by different abun¬ 
dances and different levels of approximation. Because we do not 
know the abundances in the ISM surrounding A Cephei, we use 
here the analytic representation by Siewert et al. ( 2004| l, which 
lies in between all the other cooling functions mentioned above, 
and can therefore be considered as a useful principal representa¬ 
tion. 

Mainly the hot shocked ISM is affected, the stellar wind gas 
is not. The reason for this is that the number densities in the 
shocked ISM are higher by a factor of 3.6, which is the com¬ 
pression ratio between the shocked and unshocked ISM, and that 
it is by a factor 10^ hotter, that is, TisM.shocked ~ 2.5 ■ 10^ K. 
In this region most of the mentioned cooling functions yield 
similar values. Inside the inner astrosheath, that is, the region 


M [Mo/yr] 

1.5- 10-'’ 


terminal velocity [km/s] 

2500 


spectral type 

06If(n)p 


distance [pc] 

649+^2 


radial velocity (redshift) [km/s] 

-75.1 


parallax [mas/yr] 

-7.46 and 

-11.09 

1. Parameters for A Cephei |van Leeuwen 

2007 1 . 



at 0.03 pc ISM 

n [part./cm ’’] 

6 11 

V [km/s] 

2500 80 

T[K] 

917 9000 

Rts [pc] 

0.66 

Rap [pc] 

0.93 

Rbs [pc] 

1.1 


Table 2. Stellar and interstellar boundary conditions. The stellar¬ 
centric distances of the termination shock Rts, the astropause 
Rap, and the bow shock Rps are from the model described here. 


The stellar wind parameters are taken from Henrichs & Sudnik 


(2013 1 . As interstellar temperature we use the diffuse ionized 
medium temperature derived by Reynolds et al.| ( [T99^ . The den¬ 
sity value is estimated from integrating the HI slice at the postion 
and the velocities expected for the distance of A Cephei based on 
data of the Canadian Galactic Plane Survey (Taylor et al. 2003|l. 


between the AP and the TS, the number density is on the or¬ 
der of n = 10“^ cm“^ and the velocities are on the order of 

V = 600km/s, and thus the cooling length scale, depending on 

V and n, becomes huge and is neither important there nor inside 
the TS, see below. 

The heating length scale only depends on the number den¬ 
sity n and also increases to scales much larger than the distances 
inside the AP. Therefore, as a result of the huge ram pressure of 
the stellar wind inside the TS, we can, from a dynamical point 
of view, safely neglect heating and cooling, because it does not 
influence the adiabatic expansion of the stellar wind. Moreover, 
the thermal pressure is negligible compared with the ram pres¬ 
sure in the inner astrosphere. With the help of the momentum 
equation, we can uniquely determine the shocked thermal pres¬ 
sure, which dominates in the inner astrosheath. 

In the following we always assume quasi-neutrality tig = 
tip = ti, where tig, tip are the electron and proton number density, 
respectively. Furthermore, we neglect the contribution of heavy 
ions. In Table we summarize the parameters of A Cephei de¬ 
rived from observations and the characteristic distances of its 
model astrosphere. The parallax translates into a tangential ve¬ 
locity of 41.1 km/s, which together with the radial velocity of 
75.1 km/s gives a total speed of 85.5 km/s. Our hydrodynamic 
model is three-dimensional (for later use), and in view of the un¬ 
certainties of the ISM state, we have chosen a set of parameters 
as given in Table The standard procedure in modeling is to 
choose one axis as the inflow direction (here the y axis). As long 
as the ISM is homogeneous, with a vanishing magnetic held, 
and the stellar wind is spherically symmetric, the astrosphere is 
symmetric around the inflow direction. 

The derived inner boundary conditions for the model are es¬ 
timated from the mass-loss rate and the terminal velocity and 
taken at 0.03 pc. They are presented in Table 3 together with 
those of the ISM. 

As stated above, the cooling acts differently in the subsonic 
regions from the way it does in the supersonic regions because 
in the former the thermal pressure P is dominant, while in the 
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r[K] 

V [km/s] 

n [#/cm 

^cool.s 

“^cooUs 

^heat,s 

'^heat,s 

ISMj: 

K? 

20 

40 

25kAU~ 0.12 pc 

6kyr 

248 AU « 10-3 pc 

58 kyr 

SWj: 

10 ® 

10^ 

10-3 

51GAU^250kpc 

242 Myr 

1241GAU~6Mpc 

5.8 Gyr 





^cool.h 

'^coolji 

^heat,h 

^ heat,h 

ISMi 

10 ^ 

80 

10 

56kAU^0.3 pc 

3.3 kyr 

2.8 MAU^ 14 pc 

168 kyr 


Table 3. Some characteristic numbers using the cooling function from 


Siewert et al. 


(20041. 


latter it is the ram pressure pv^ 12, where p and v are the mass 
density (mainly protons, but helium or other elements may con¬ 
tribute to the mass density) and bulk speed, respectively. While 
in the following we discuss only protons, we can add, for ex¬ 
ample, helium, which leads to partial densities, temperatures, 
and pressures for which the approximations made below can be 
separated. We are interested in the shortest characteristic length, 
which is in the subsonic case inversely proportional to the num¬ 
ber densities, and thus the protons dominate. 

In the hypersonic case, the estimation below can differ by up 
to 40% when including helium. This would then also require in¬ 
cluding helium as a new species in the Euler energy equations 
(including different species, see [Scherer et al. |2014| l, which 
would violate our assumption of a single fluid. The interaction 
of other species concerning the energy loss by Coulomb colli¬ 
sions is, however, already included in the cooling functions, so 
that we can continue with the single-fluid equations consisting 
of protons including the cooling term for a first analysis. 

From the stationary energy equation we obtain in the sub¬ 
sonic region with the assumption nmpi?l2 <Si P: 


V- 


y 1 

- -P + -nm„v \v 

y - 1 2 


jPv 


iy ^)^cooi,s 


( 1 ) 


5nkT V 

^cool,s 


= -n^MT) 


where y = 5/3 is the adiabatic index, irip the proton mass, 
P - 2nkT, k is the Boltzmann constant, v, T the bulk velocity 
and temperature of the plasma flow, and A the cooling func¬ 
tion. Taking the absolute values in Eq. Q and replacing V by 
the inverse subsonic cooling length, Lcooi.s can be estimated as 
follows: 


^cool,s 


yPv 


5kTv 


(y-l)n^A{T) nA(7’)’ 

and the subsonic cooling time TcooI/- 
^cooLs 5kP 


^cool,s 


nK{T) ■ 


( 2 ) 


(3) 


This cooling time is up to a factor 3 the same as that given in 
Sutherland & DoprE| ( |1993l l. 

For supersonic interstellar flows, that is, PismV^^^/2 » y/(y- 
l)P, it follows: 


y 1 

—P + -niUpV" 


v^O 


MptlV 

2l 


'cooLh 


(4) 


and we derive the supersonic (index h for “hypersonic” to dis¬ 
tinguish it from the subsonic one) cooling length Lcooiji and time 

^ cooLh 


^cool,h 


2nA(T) 


and 


'^cooLh — 


2nA(T) 


(5) 


For the heating function by photo-ionization and some sup¬ 
plemental heat source we use the approach by [Reynolds et al.j 


(1999 1 ; see also Kosihski & Hanasz ( 2006| l. The heating rate for 
photo-ionization depending on electron collisions is limited by 
recombination and is thus proportional to n^, while additional 
heating terms that can include photoelectric heating by dust, dis¬ 
sipation of turbulence, interactions with cosmic rays are propor¬ 
tional to 


F = nGo + nG\, 


-25 , 


10 ^"^erg cm^ 


( 6 ) 

' and Gi = 


with the co nstants Gn - 

' (Kosihski & Hanasz[ |2006|). Replacing the right 


10 erg s 
hand side of Eq. ([IJ by F, we obtain 


^heat,s ~ 


5kT V 
nGo + Gi 


and 


Lh, 


eat,h 


fripV 


2(nGo + Gi) 


(7) 


and the heating times Theat,s,Theat,h by dividing the respective 
heating lengths by the speed v. 

We can estimate the cooling and heating lengths and times 
for the shocked ISM (ISM 2 ), and analogously for the shocked 
stellar wind SW 2 (see Fig. [^. The results are displayed in 
Table together with the characteristic lengths and times for 
the hypersonic ISM (ISMi), because it is also cooled. For the 
ISM the dependence of the cooling and heating lengths is shown 
in Fig. 1^ The vertical line denotes the temperature at which the 
ram pressure equals the thermal pressure. Feft of this line the hy¬ 
personic length scales from Eqs. (j^ and (0 are displayed, while 
on right the subsonic scales are shown (Eqs. (|5 and (0). 

The supersonic parameters are number density n = 11 cm~^ 
and a speed of y = 80 km/s while for the subsonic parameters, the 
density was multiplied by the compression ratio i = 3.62 and the 
speed was divided by s using the Rankine-Hugoniot relations. 

Figure shows that in the subsonic case for temperatures 
above 5 ■ 10^ K the heating scale length is always longer than 
that for the cooling. Thus cooling is more efficient tfian fieat- 
ing for the discussed functions and parameters. In the supersonic 
case (T < 5 ■ 10^ K), the heating scale lengths are only longer 
down to temperatures of =» 10"^ K, while for temperatures below 
^ 4 ■ 10^ K the cooling scale length is more important. Thus, 
above ^ 4 ■ 10^ K cooling is more efficient tfian heating, and 
below this, heating is important. We can read from this figure 
tfiat the length scale for cooling of the shocked ISM tempera¬ 
ture T ss 2.5 ■ 10^ K is ss 0.01 pc. From the model we see that 
the BS to AP distance is 0.17 pc, which is more than ten times 
the cooling length Lcoo/,i- This strong cooling is balanced by the 
fact that at least at the stagnation fine tfie sfiocked ram pressure 
l/ 2 p 2 ,/SA/t^ 2 «/SM converted into tfiermal pressure to¬ 

ward tfie astropause, througfi wfiicfi mass flux is zero. 

It is also evident from Fig. [^that if the temperature falls be¬ 
low ^ 4 ■ 10^ K, the flow again becomes supersonic. 

Different values of the number density n or speed v or dif¬ 
ferent cooling or heating functions change the characteristic 
lengths, but the principle estimates for the dynamics of the flow 
field remain the same. 

In the stationary case, where no relative speed between the 
star and the ISM exists, one should take for the speed v that of 


3 










































K. Scherer et al.: Cosmic rays in astrospheres 



2.0 3.0 4.0 S.O 6.0 7.0 8.0 9.0 

loglO (T[K]) 


Fig. 2. Number density and speed are fixed at n = 11 part/cm^ 
and V - 80km/s. All scales are logarithmic. On the horizon¬ 
tal axis the temperature is displayed, while on the left vertical 
axis the characteristic length is shown in AU and in parsec on 
the right y-axis. The black vertical line denotes the temperature 
where the ram pressure equals the thermal pressure. 


the outward-moving shock front. Moreover, the criterion given 
by Schwarz et al. (19721, namely that £/log(A(7’))/£/log(7’) > 2 
to guarantee thermal stability, is generally not satisfied, and thus 
one expects clumping in between the AP and the BS and possibly 
in the ISM. The cooling lengths can be scaled by the speed and 
the density: Lcooi,s does not change as long as u/n = const. This 
is not true for the heating length Lheat,s- 


A similar result was also derived by van Marie & Keppens 
( 201 l| l. In the supersonic ISM the ram pressure dominates the 
cooling scales. It is expected that a magnetic field will increase 
the characteristic length scale due to its additional pressure and 


thus helps in stabilizing the shock structures ( van Marie et al 
2014|l. The above estimate of the characteristic cooling length 


is crucial in determining the resolution required for numerical 
models of the large-scale structure of astrospheres. From Fig.[^ 
it is clear that this should be higher than the shortest charac¬ 
teristic cooling lengths, which are on the order of >0.01 pc for 
the numbers given there. Moreover, the cooling of the surround¬ 
ing ISM can also be easily be inferred from the above estimate, 
and thus it can be approximately determined whether it is cooled 
during the computation. 


4. Astrosphere model 


The hydrodynamic model for the star A Cephei is made in 3D; 
the third dimension is needed for future comparison with models 
including bipolar winds or magnetic fields. The boundary con¬ 
ditions are given in Table While the stellar mass loss and the 
terminal speed can be determined by observations, the stellar 
wind temperature as well as the interstellar parameters are so¬ 
phisticated guesses. The model solves the Euler equations for A 
Cephei using the Cronos MHD code as described in Kissmann 
et al. ( |2008 1 , Kleimann et al.| ( |2009) l, and |Wiengarten et al. 
( 2013| l. The results are displayed in Fig.j^for the proton number 
density of A Cephei. In Fig.j^the wiggles along the AP caused 
by the thermal instabilities can be clearly recognized. They are 
due to the cooling functions. 

This model provides the underlying plasma structure needed 
as input for the transport equation discussed below. The model is 
solved on a spherical grid with a resolution of 0.005 pc in radial 
dimension, and 2° and 3° in § and ip dimension, respectively. The 


bg( n ) [cm''-3] 



-10 -50 5 10 

X [pc] 


Fig. 3. Logarithmic number density for A Cephei (including 
heating and cooling). The axes are given in pc, while the color 
bar is a logarithmic scale for the proton number density. 



Fig. 4. Continuum-corrected and slightly spatially filtered Ha 
image of the bow shock around A Cephei. North is up, east is 
left, the size of the plotted field is 22' x 22'. The imaging data 
were taken as part of the IPHAS ( |Drew et'af] 2005| l survey. A 
bright and well-defined bow shock nebula is visible toward the 
SE of the star. The roughly spherical Ha structure centered on A 
Cephei is an out-of-focus ghost of the bright star, the dark inner 
rings are the equivalent structure from the r-filter image used for 
continuum subtraction. In addition to these purely instrumental 
effects, several Ha bright features are visible close to the star. 
A more detailed discussion of the multiwavelength properties of 
the bow shock and the circumstellar environment of this runaway 
06 star will be presented in an upcoming paper. 


large distance of the outer boundary is needed because during 
the evolution of the astrosphere it becomes much broader than 
shown in Eig.j^and finally shrinks to the state shown here after 
ca. 170kyr. 
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5. Cosmic ray fluxes 

To model the flux of OCRs, the Parker transport equation 


If =-i'-V/ + V-(K-V/) + 


df 


3 dlnP 


has been frequently used in the literature (see Potgieter et al. 


2001 and reference therein). In this equation, / is the (nearly) 


isotropic GCR distribution function, v the bulk plasma flow, K 
the diffusion tensor (including, in a 3D geometry, separate com¬ 
ponents directed parallel and perpendicular to the mean mag¬ 
netic field), and P is the particle rigidity. As an initial approach, 
spherical symmetry is assumed, so that the Parker equation re¬ 
duces to 


dfirj) df 1 d 


dt 


dr 


2 ^ 
r^ dr \ dr 


P d , 

3r2 dr 


df 

51nP’ 


(9) 


where r is radial distance and Krr is the effective radial diffusion 
coefficient. In this work, Eq. (j^ is solved by transforming it into 
the equivalent set of stochastic differential equations (SDKs) 


dr — 


dP = 


1 d 


^ ( 2 \ 
r2 /ir V 


'■ dr ^ 
P d 
3r2 dr 


dt + yjlKrrdW, 


('‘o) 


dt. 


( 10 ) 


( 11 ) 


which is then integrated numerically (Strauss et al. 201 l| l, where 
the solution of the ID equation along r can be found in|Strauss| 
[etar] ( |20TT] l. 

These equations are coupled to the simulated HD geometry 
by reading in the modeled values of v and V -v (governing energy 
changes) directly from the HD simulations (along the stagnation 
line) and solving for the GCR flux using essentially a test par¬ 
ticle approach. The magnetic field enters the computations via 
the diffusion coefficient Krr - Krr(B). In the ID scenario, the 
geometry of B does not enter the computations, although for a 
azimuthal field (the case inside the TS), Krr ~ thus reducing 
to a diffusion coefficient perpendicular to the mean field. As a 
boundary condition for the GCR flux, a local interstellar spec¬ 
trum is specified at the edge of the computational domain. The 
GCR differential intensity is related to the distribution function 

by i = P^f- 

Based on experience gained from modulation studies inside 
the heliosphere, Krr can be decomposed into a radial and energy 
dependence, Krr - {r)K 2 iP), where, for this study 


Kiir) = 


Ksw if ?• < rss 
^ism if r > rj2s 


( 12 ) 


with rss the radial position of the BS and where Kjsm is inde¬ 
pendent of position and Ksw = kqBqIB. kq is a normalization 
constant, usually specified near the inner boundary, and Bq is a 
constant added for dimensional consistency. Assuming that the 
astrospheric magn etic field is about 80 times higher than the 


heliospheric case (Naze 2013 i, a value of 0.027 kpc^Myr ‘ is 


used, scaling up the heliospheric diffusion coefficient by about 
two orders of magnitude for this astrospheric case. The B ' ra¬ 


dial dependence is based on the results of Engelbrecht & Burger 
( 2013| l, but approximated in such a fashion because in situ ob¬ 
servations are, of course, not available. The energy dependence 
of Krr is taken from Biisching & Potgieter (2008|l 


KiiP) = 


+ 0.6 


— if P> P, 


0 


0.48 


if P< Pi 


0 


with Po - 4 GV. Instead of showing Krr, we follow the conven¬ 
tion of showing the corresponding mean free path (mfp). 


3Krr 

( 8 ) , 


(14) 


where Vp is particle speed. The magnitude of Krr is changed in 
the next section to illustrate its effect on the resulting particle 
intensities. 

Because the mfp used is, in general, not known, we studied 
the behavior of three different values inside A Cephei and three 
different values of the interstellar mfp. 

5. 1. Mean free path 

The structure of the astrosphere can be recognized in the lower 
part of Eig.j^by the respective jumps in the mfps (from the right, 
which is the inflow direction): the BS (1.1 pc), AP (0.93 pc) and 
TS (0.66 pc), represented by the vertical dotted lines. This TS is 
marked by sharp decrease of the mfp, the AP by a change of the 
slope, and finally the BS by the sharp drop at the ISM side. The 
resulting differential flux (DI) for 1 GeV particles along the stag¬ 
nation line is displayed in Pig.|^ In all cases the modulation, that 
is, the decrease in the differential intensity DI, starts far away in 
front of the BS, a feature that was discussed by IScherer et al. 


( |201 l| l, [Herbst et al.| ( |2012| l, and [Strauss et al.| ( 2013| ) for the 
case of the heliosphere. 

This outer astrospheric modulation depends on the ratio be¬ 
tween the parallel and perpendicular diffusion coefficient and 
vanishes for high ratios. Thus, for an appropriate choice of the 
transport parameters, the modulation of GCRs starts in front of 
the astrosphere and then rapidly decreases at the astrospheric BS. 
It is also evident that in contrast to the heliosphere, the GCRs are 
modulated directly behind the BS: This is due to the cooling ef¬ 
fects, which shrink the region between the AP and BS to 0.17 pc 
and are a barrier-like feature for the cosmic rays ([Potgieter &| 
[Langner|[2004[ ). 

In the lower part of Pig. we plot the mfp used in the mod¬ 
eling. The inner mfp (Asw) obtained from the model are then 
divided or multiplied by a factor two and three to demonstrate 
its influence. A short Asu, strongly increases the modulation, 
and the flux of 1 GeV particles almost vanishes inside the AP 
. Increasing Asu, by the same factor leads to a nearly vanishing 
modulation. Changing the interstellar Ai^m from 50 pc to 100 pc 
and 500 pc does not change the DI remarkably. These variations 
show the theoretical expected effects. Especially for short mfp 
inside the astrosphere the GCR spectra are efficiently cooled. 
The 20% modulation in front of the astrosphere is presumably 
caused by the sharp drop of the mfp between the AP and the BS. 

Prom Pig. 1^ we can also see that most of the modulation 
- up to a factor 5 or even more - occurs between the BS and 
the AP (depending on the transport parameters). Thus, because 
astrospheres are large volumes surrounding the star, the GCR 
spectrum can be cooled and, especially, small-scale anisotropies 
in the galactic CRP can be observed far away from large astro- 
spheres. This prediction will be improved by further modeling 
and studying the distribution of hot stars in the Galaxy. 

At the BS, the mfp falls to quite low values, rises sharply at 
the AP, and finally jumps again at the TS. In the inner astrosheath 
the mfp increase continuously toward the shock because we as¬ 
sumed that the magnetic field is inversely proportional to the 


(13) speed v as for the Parker spiral field in the heliosphere (Potgieter 


et al. 2001 1 . Thus at the stagnation line, the speed must vanish 


at the AP because this is a tangential discontinuity with no mass 
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Fig. 5. Modulation of 1 GeV particles in A Cephei along the 
stagnation line for different stellar and interstellar mfp. The mfp 
are shown in the lower panel and the colors correspond to those 
in the upper panel. The magenta, orange, and green curves in the 
upper panel are identical when normalized to their respective 
values at infinity. The dotted lines show the position of the TS, 
AP, and BS. 


flow through it. Inside the TS, the magnetic field increases with 
r in the direction to the star, and the mfp decreases. Because of 
the strong compression of the outer astrosheath, the mfp dimin¬ 
ishes to an even lower value, which then forms a barrier for the 
cosmic ray transport. 

This qualitative behavior is theoretically well understood and 
confirmed by observations in the heliosphere, except for the ad¬ 
ditional modulation in the outer astrosheath. The latter differs 
because of the cooling function. Thus, the cooling function has 
a direct influence on the cosmic ray transport. 

5.2. Energy dependence 

In Fig. 1^ we study the behavior of the GCR flux at different en¬ 
ergies. The 1 GeV particles are modulated by a factor of almost 
five, while for higher energies this factor becomes smaller, as ex¬ 
pected. But even for the highest energies of 10.. lOOTeV, a small 
modulation of a few per mille (1 per mille = 0.1%) outside of the 
BS is visible. 

In Fig.|^we plot the DIs for different energies, from 1 GeV to 
100 TeV in logarithmic steps of ten. In the upper panel the lower 
energies (form 1 to 10 GeV) are strongly modulated, while for 
the 10 GeV and 1 TeV particles the DI is only diminished by a 
few percent (middle panel). The lower panel shows the modu¬ 
lation of the 10 and 100 TeV particles, which is in the per mille 
range or even much lower for the 100 TeV particles, which are 
hardly affected. All energies are modulated ahead of the BS, and 
thus the modulation volume is larger than that of the astrosphere. 
The DI of the TeV particles can be increased by diminishing the 
mfp inside the astrosphere or can become less by increasing the 
mfp; this is not shown in Fig. The mfp is based on that in the 
heliosphere, and the modulations displayed in Fig. are based 
on it. 

5.3. Observational evidence 

Small-scale modulations of a few per mille in the TeV range 
are observed with, for example, the IceCube experiment (e.g., 
[Abbasi et al.| |2012| and references therein). In particular, their 


0 12 3 4 



Fig. 6. The modulation of particles with different energies. In the 
upper panel are the normalized differential particle intensities in 
part.fiirP's srMeV) for energies of 1 and lOGeV are shown (in 
black and blue, respectively). The y-axes are in a linear scale 
with different values for all three panels. In the middle panel the 
DI for particles at 100 GeV and 1 TeV are displayed (black and 
blue). In the lower panel the DIs for 10 and 100 TeV particles are 
shown. All Dis are normalized to their corresponding interstellar 
DI and are decrease towards the star. 


Fig. 5 shows a resolution of 3° , which cannot resolve the astro- 
spheres of A Cephei with angular diameters 0.2". Nevertheless, 
the pixels show variations that may be attributed to local GCR 
fluctuations, provided the statistics in the pixels is large enough 
and the astrospheres are close enough. Because runaway O/B 
stars or OB associations are quite frequent in the Galaxy 
( Huthoff & Kaper] 2002| l, their local disturbance of the GCR 
spectrum can explain small-scale variations in the GCR flux. A 
Cephei was only used as an example, but the above holds true in 
general for all large-scale astrospheres. 

Our simulations may help to understand these variations. The 
variation in the Dis can be increased or decreased by choosing 
a smaller or larger mfp inside the astrosphere, therefore a better 
estimate of their magnetic fields is needed, which can help to 
understand their turbulence levels. 

For the higher energies a modulation outside the astrosphere 
can also be observed, where the “astrosphere of influence” is 
roughly twice the BS distance for bullet-shaped astrospheres like 
that of A Cephei. This holds true in the direction of the inflow¬ 
ing ISM, while it is more complicated in other directions. These 
details will be explored in future. We studied here only protons, 
but because the transport coefficients depend on rigidity, the be¬ 
havior of other species can roughly be estimated. For helium, for 
instance, a similar behavior is expected as for the protons shifted 
by a factor of two in Fig. At these high energies the rigidity 
has roughly the same value as the particle energy. This means 
that multiplying the rigidity by a factor of 2 is the same as do¬ 
ing it in the particle energy. Therefore, helium is slightly less 
modulated when passing through astrospheres than the protons. 

If there were a few large astrospheres (or stellar wind bub¬ 
bles) in the direction toward an observer, the GCR spectrum 
would be slightly cooler than in other directions. This can 
contribute to the understanding of the small-scale anisotropies 
present in the IceCube data (Abbasi et al. 2012| l. 

The parameter set for the supersonic wind and ISM are cho¬ 
sen to be close to the observed values, but for modeling purposes 
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we worked with rounded values. In addition, for the transport 
model we took the turbulence levels based on that of the he¬ 
liosphere. Both of these aspects need further attention, but the 
principal effect remains: large stellar astrospheres can affect the 
local interstellar cosmic ray spectra. 

6. Conclusion 

Based on our models, we studied the transport of OCRs in an 
astrosphere. We have shown that even ahead of an astrosphere, 
where the mfp is still undisturbed, a cooling of the GCR spectra 
occurs because particles are trapped by scattering into the astro- 
sphere, in which they can be effectively cooled. The “astrosphere 
of influence” around the modeled astrospheres is roughly twice 
its hydrodynamic dimension. Thus, stellar wind bubbles (astro- 
spheres) can cool the Galactic spectrum, and small anisotropies 
are expected in the direction to an observer. 

As a result of its effect in compressing the outer astrosheath, 
the cooling function directly influences the GCR modulation in 
this region. Additionally, the cooling and heating functions re¬ 
quire high resolutions for global models because of their char¬ 
acteristic length scales. They also give a rough estimate whether 
the surrounding ISM is also affected during the calculations. 

The modulation in astrospheres affects particles up to 
lOOTeV. This can help to understand the anisotropies on small 
angular scales, for example, by the IceCube experiment, among 
others. The angular extent of A Cephei is too small to be resolved 
by these experiments because of its large distance (ss650pc). 
Nevertheless, nearby hot stars or stellar associations can have 
an larger angular extent and may possibly be observed. In our 
simulations we saw an effective cooling of the GCR spectrum. 
Thus, the conclusion is that small-scale cosmic ray anisotropies 
may be explained by the modulation in such huge cavities. 

Because the model is fully 3D, the modulation or other pa¬ 
rameters along a line of sight toward Earth can in principle be 
theoretically determined. Here we demonstrated that large as¬ 
trospheres can modulate the Galactic cosmic ray spectrum quite 
significantly, and no homogeneous spectrum over all directions 
can be expected. To calculate the modulation of GCRs along a 
line of sight or at Earth requires more sophisticated methods than 
discussed here, but this is being developed. 
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